library("sf")
library("raster")
library("rasterVis")
library("ncdf4")
library("readxl")

# perform data diagnostics 

sst <- stack("input/cmems_mod_blk_phy-tem_my_2.5km_P1M-m_1715770301759.nc")

### Visualize anual SST by all dataset years
at <- c(10, 12, 14, 16, 18, 20)
cc <- list(at=at, ## where the colors change
           labels=list(
             at=at## where to print labels
           ))

years <- 1993:2021
for (y in years) {
  i <- y-1993
  
  offset <- (i*12) + (1:12)

  slice <- sst[[offset]]

  png(paste0("output/mapsst/", y ,".png"), res=180, width=1200, height=900)
  print({
    levelplot(mean(slice),  contour = F, region = T, margin = F, main = paste0("Среднегодовая ТПВ ", y), par.settings = magmaTheme, at=at, colorkey = cc)
  })
  dev.off()
}

### perform ACF & boxplot
data <- read_excel("input/dataset_prepared.xlsx")

# anchovy
acf(na.omit(data$anchovy_nw_catch), lag.max = 5)
boxplot(na.omit(data$anchovy_nw_catch))
shapiro.test(na.omit(data$anchovy_nw_catch))

# bonito
acf(na.omit(data$bonito_catch), lag.max = 5)
boxplot(na.omit(data$bonito_catch))
shapiro.test(na.omit(data$bonito_catch))

# anual sst
acf(na.omit(data$sst_anual), lag.max = 5)
boxplot(na.omit(data$sst_anual))
shapiro.test(na.omit(data$sst_anual))

# winter sst
acf(na.omit(data$sst_winter), lag.max = 5)
boxplot(na.omit(data$sst_winter))
shapiro.test(na.omit(data$sst_winter))
